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We present a new, computationally inexpensive method for the calculation of reduced density matrix dynamics for 
systems with a potentially large number of subsystem degrees of freedom coupled to a generic bath. The approach 
consists of propagation of weak-coupling Redfield-like equations for the high frequency bath degrees of freedom only, 
while the low frequency bath modes are dynamically arrested but statistically sampled. We examine the improvements 
afforded by this approximation by comparing with exact results for the spin-boson model over a wide range of parame¬ 
ter space. The results from the method are found to dramatically improve Redfield dynamics in highly non-Markovian 
regimes, at a similar computational cost. Relaxation of the the mode-freezing approximation via classical (Ehrenfest) 
evolution of the low frequency modes results in a dynamical hybrid method. We find that this Redfield-based dynam¬ 
ical hybrid approach, which is computationally more expensive than bare Redfield dynamics, yields only a marginal 
improvement over the simpler approximation of complete mode arrest. 


I. INTRODUCTION 

Useful approximate methods for the description of quan¬ 
tum dynamics and relaxation can often circumvent the large 
computational expense of numerically exact approaches while 
maintaining quantitative accuracy in certain regions of param¬ 
eter space. The general applicability of such methods, how¬ 
ever, is often limited and their domain of validity difficult to 
assess. The most widely used approximate approaches fall 
into two broad and general classes. The first class of methods 
involves techniques that employ uncontrolled approximations 
to yield dynamics which are non-perturbative in the various 
couplings (e.g., intra-system or system-bath) that character¬ 
ize the problem. The second class of methods are systemati¬ 
cally perturbative in a well-defined coupling parameter, but 
are free from further classical or semiclassical approxima¬ 
tions, at least for simply defined models such as a the spin- 
boson model. 

One of the most celebrated perturbative techniques is a 
lowest-order treatment of the system-bath coupling, known 
traditionally as Redfield theory. 13 As we will show later, the 
relevant dimensionless parameter characterizing the accuracy 
of Redfield theory is r/ = max 1 2A/flotj, 2A/na> c |, where A is 
the nuclear reorganization energy, fl = 1 /k B T is the inverse 
temperature, and u> c is a characteristic bath frequency (we 
henceforth work in units with ft = 1). Redfield theory be¬ 
comes unreliable when t) > 1. We emphasize that // is con¬ 
trolled by multiple bath parameters and, in particular, low fre¬ 
quency degrees of freedom (small a> c ) limit the range of acces¬ 
sible reorganization energies. Indeed, violations of this condi¬ 
tion explain the failures of Redfield theory found by Ishizaki 
and Fleming for certain models of excitation energy transfer, 
which appear to be characterized by low-frequency protein 
baths. 

While the very lowest frequency degrees of freedom are 
thus most problematic for Redfield theory to handle (even 
in its non-Markovian forms), it is often the case that nu¬ 
clear modes of such frequencies are effectively frozen on the 


time scale of relevance for the system’s dynamics. In this 
regard, the key function of such modes is simply to pro¬ 
vide static energetic disorder for the more rapidly evolving 
degrees of freedom. This suggests a methodology whereby 
the very low frequency phonons are approximated as static 
(and treated non-perturbatively as a source of static disorder), 
while the remaining portion of the bath is treated dynamically 
within Redfield theory. Here we develop this “Redfield theory 
with frozen modes” (Redfield-FM) method, and show that it 
greatly extends the applicability of Redfield theory into highly 
non-Markovian dynamical regimes at essentially no change in 
computational cost. 

The outline of this paper is as follows. In Sec. [II] we in¬ 
troduce the theoretical background for the Redfield equations 
and the derivation and general properties of the Redfield- 
FM extension. In addition, we also introduce the spin-boson 
Hamiltonian as the model system on which we test the meth¬ 
ods developed in this paper. Section III A| presents the com¬ 
putational details in the implementation of the Redfield-FM 
method, while Sec. III B| presents illustrative results for the 
method. In Sec. IIIC we relax the mode-freezing approx¬ 
imation via the derivation and implementation of a dynami¬ 
cal hybrid method (hybrid-Redfield) that combines Redfield 
dynamics for the high-frequency part of the bath coupled to 
the electronic system and Ehrenfest dynamics for the low- 
frequency modes. In Sec.|IV| we conclude. 


II. THEORY 
A. Model 

First, we briefly describe the model system we use to test 
the approximations developed in subsequent sections. This 
allows us to define notation and parameters that will be used 
in our numerical comparisons. We focus on the well-known 
spin-boson model, which consists of a two-level system cou¬ 
pled linearly to a harmonic bath. This model has been exten- 
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sively used to investigate a wide variety of relaxation, charge 
and energy transport processes in condensed phase systems. 5 

The total Hamiltonian is divided into system, bath, and in¬ 
teraction components, H = H sys + // ha th + V. The system 
Hamiltonian takes the form 

H sys = scr z + Acr x , (1) 


e Oe' ( ' Hs ^ +Hid ^', and the free bath correlation func¬ 

tion is given by 


C (0 = ^ C l Tr bath {pbathQ*(0Q*(0)} 
k 


1 

n 



dcoJ(u))[coth(j3u>/2) cos(cot) - i sin(tuf)]. 


( 6 ) 


where cr ; , i = [x,y, z], are the Pauli matrices, 2s is the energy 
difference, and A is the coupling between the two electronic 
sites, which is here assumed to be static. The bath Hamilto¬ 
nian consists of an infinite set of harmonic oscillators, 

( 2 ) 

k 

Lastly, the system-bath interaction couples the electronic 
states linearly to coordinates of the bath oscillators, 

V-rz^cuOk. (3) 

k 

Physically, the system-bath coupling acts as a (quantum) fluc¬ 
tuating field that shifts the origin of the bath harmonic oscil¬ 
lators by a magnitude that depends on the system’s electronic 
state and the strength of the coupling. 

The spectral density, which completely determines the cou¬ 
pling between the bath and the system, is taken to be Ohmic 
with a Lorentzian cutoff (Debye form), 

c 2 

J(w) = g V —S(a> - 0 J k ) = 2Aeo c ^ (4) 

2 X-U (Ok cl>- + coi 

k L 

The cutoff frequency, co c , characterizes how quickly the bath 
relaxes toward equilibrium, while the reorganization energy, 
A = 7T~ l J dco J(cj)/co, characterizes the energy dissipated 
by the environment after a Franck-Condon transition between 
electronic states. It is important to note that the methods stud¬ 
ied here are neither limited to the spin-boson model nor to the 
Debye form for the spectral density. 


B. Time-local Redfield Dynamics 

Because of its simplicity in the time domain we employ the 
time-local (i.e. time-convolutionless) form of the generalized 
Redfield equations. A full derivation of these equations is con¬ 
tained in Appendix [A] Here our aim is to highlight important 
but often overlooked aspects pertaining to the applicability of 
the Redfield approach. For the spin-boson model, the time- 
local version of the Redfield theory takes the following form, 

-^PO) = -i[#sys.p(0] 

dr {C(t) [cr z (0), (T z (-T)p(t)\ ( 5 ) 

- C*(t) [it z (0), p{t)cr -(-t)] }, 

where all operators except the reduced density matrix 
(RDM) are evolved in the interaction picture, O(t) = 



By going to the interaction picture with respect to H sys + //bath 
(to eliminate the free-evolution) and formally integrating the 
equation of motion |(5j, one finds 


Pi(t) = P/(0) 



dT 2 C(r 2 ) 



( 7 ) 


Examination of the function f 0 'dr i fj 1 dr 2 C(t 2 ) reveals 
the natural dimensionless parameter that determines the 
limit of validity of Redfield theory. In general, even in 
the non-Markovian case, we expect that the ordering of 
terms in the expansion 0 is governed by the function 
f () dr\ fj 1 dr 2 C(r 2 ) = r]g(t), where 7/ is a dimensionless con¬ 
stant and g(t) is a function expressed in terms of a scaled, 
dimensionless time variable. In the high-temperature limit 
( J5io c <s 1), where C(t ) * (2A//3)e^ 0>c ', it is easy to show 

i] = 2 A/(/3 co 2 c ) (8a) 

gif) = e- 1 + cop. (8b) 


In the low-temperature limit iJ3u> c » 1), we assume that the 
low-frequency behavior of the spectral function dominates, 
so we choose J(a>) = (2 Aco/co c )e^ 2a, ^ 0,c as an approximation 
to the Debye form in Eq. (|4]> that exactly matches the value 
of /l and its low-frequency asymptotic behavior. Using this 
spectral density, one can show 



ncj c 


git) = ln[l + (to c tfj + 2 In 


smhint Ipti) 
(nt/pti) 


(9a) 

(9b) 


Thus, for a Debye spectral density, Redfield theory will be 
reliable as long as max ^2A//3(op2A/mOc^ is not significantly 
larger than unity.® It should be noted that recent work pur¬ 
ported to be in the Redfield limit actually violates the above 
condition. 7 As long as the relevant energy scales in the system 
Hamiltonian are not too large, we expect the above to hold. In 
cases where the system’s bare energy difference s is the largest 
energy scale in the problem, the dynamics will be mediated 
by multi-phonon processes which are a challenge for lowest- 
order Redfield-like theories. However, in this limit, the prob¬ 
lem acquires an increasing amount of ‘pure-dephasing’ char¬ 
acter, for which the time-local version of Redfield theory pro¬ 
vides an exact multi-phonon resummation. 


C. Redfield Theory with Frozen Modes 

As discussed in the Introduction, low-frequency bath 
modes co^ which lead to a violation of the validity of the Red¬ 
field theory frequently evolve so slowly as to be effectively 
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static on the electronic timescale. Here we develop the “Red- 
field theory with frozen modes” (Redfield-FM) method, based 
on the physically appealing notion of dividing modes into a 
low-frequency portion (treated as static disorder), and a high- 
frequency bath (treated by time-local Redfield theory). The 
separation of modes used here mirrors that utilized in previ¬ 
ous work on a Forster-like dynamical hybrid approach.® The 
approach presented here is not in any way limited to time- 
local Redfield theory (nor even to any specific flavor of Red- 
field theory). However, due to pathologies associated with a 
strictly Markovian Redfield theory, we suggest certain adjust¬ 
ments to the partitioning algorithm, as discussed in Appendix 

m 

First, it is advantageous to assume that the total density 
matrix is multiplicatively separable into weakly interacting 
parts, i.e., p tot (f) * p s iow(0)Psy S +fast(f) where p s i ow (0) is the 
density matrix of the frozen low-frequency “slow modes” 
and p S ys+fast(f) is the density matrix for the system and high- 
frequency “fast modes”. As in Ref. [8j a splitting func¬ 
tion, S (co), divides the spectral density into two components, 
J(oj) = /siow(w) + Tfast(w), where 

J s1ovi (oj) = S(oj,oj*)J(oj), (10) 


and 


(a) of = 0.9cj c 


(b) a>* = 2.0u> c 



FIG. 1. Spectral density and splitting via S {cj, w”), illustrating the two situa¬ 
tions expressed in Eqs. 0 and 0 


Vslow = CT Z Efeslow CkQk ■ Regrouping terms, it is evident that 
freezing the slow part, V s i ow , will yield a classical reorganiza¬ 
tion energy that renormalizes the bias for every realization of 
the bath’s initial conditions. The modified total Hamiltonian 
now takes the following form, 

H' = [s + A cl (0)]cr z + A cr x 

+ a z^ c k Qk+2^[ P k + (0 kOk\ <13) 

kefast k 


•W«>) = [1 -S(0J,0/)]J(0J) (11) 

Here we take the same form of the splitting function as that 
suggested in Ref. [8] namely. 


S{u,of) 


[1 - ( tn / w *) 2 ] 2 

0 


oj < of 
o> > of, 


( 12 ) 


which, by virtue of its smoothness, avoids problems associ¬ 
ated with long-time oscillatory tails in the bath correlation 
function.® While the above splitting induces no errors if the 
dynamics are treated exactly, it is clear that of serves as a free 
parameter that allows one to tune the optimum percentage of 
frozen bath modes, and hence the accuracy of the results if 
the dynamics are treated within our approximate method. The 
utility of the present method is greatly enhanced if a physical 
a priori prescription for choosing of based only on the pa¬ 
rameters of the initial Hamiltonian can be put forth. In this 
work we choose of = ojr/4, where ojr — 2 ye 2 + A 2 is the 
system Rabi frequency. This choice for of is simple, yields 
non-trivial improvements over standard Redfield theory, and 
may easily be generalized to multiple electronic states. Phys¬ 
ically, this choice partitions the bath into modes that evolve 
slower than the system (to be treated as frozen) and modes 
that evolve faster than the system (to be treated via Redfield 
theory). However, it should be noted that this choice is not 
always optimal. Future work will be devoted to the goal of 
arriving at an optimal choice of of. Other choices for of that 
fit within the general physical guidelines discussed above will 
be discussed in the Results section. 

With a prescription for choosing of in hand, it is possi¬ 
ble to separate the fast and slow portions of the Hamilto¬ 
nian, starting with the interaction, Vf ast = cr, £ i6fast c kQk and 


where the classical reorganization energy is defined as 
T c; (0) = Xligsiow CkQk(0) and the set of Qk(0) is sampled from 
a bath distribution function after the discretization of ./ s | ow (w). 
Physically, each realization of the frozen bath degrees of free¬ 
dom constitutes a local, rigid environment that modifies the 
site energies for the system Hamiltonian. The time-local Red¬ 
field dynamics, under the Hamiltonian in Eq. ( [T3] l with an in¬ 
teraction given by Vf ast = cr. 2 fefast cuQk, are subsequently en¬ 
semble averaged with respect to the slow frozen modes. Thus, 
there are two important differences for the Redfield equations 
used in each realization of the frozen modes: (i) the bias is 
given by s = s + A cl ( 0), and (ii) the bath correlation function 
given by Eq ([6]l is modified, with J(oj) replaced by 7f ast (w). 

To account for the classical frozen modes, the nonequilib¬ 
rium population dynamics takes the following form, 

(cr,(f)) ~ I dPdQ p s iow(P, Q, 0) Tr sys+ f as t[psys+f a st(0<T,], 

(14) 

where p s i ow (P. Q, 0) could be either the classical distribution 
function or the Wigner transform of the equilibrium density 
operator of the slow bath degrees of freedom, and p S ys+fast(0 is 
the reduced density matrix of the system and the fast bath de¬ 
grees of freedom. Observables such as Tr sys+ f as t[p s y s+ f as t(f)cr;,] 
may then be calculated via the Redfield equations. 

To understand the relaxation processes with mode freezing 
for finite of , it is first useful to investigate the effect of the 
approximation at its most extreme, namely the adiabatic limit, 
where all bath modes are assumed to be static {of —» oo). 
In this limit, we are effectively in the Born-Oppenheimer 
regime, where excitations in the electronic subspace move 
along the potential energy surface determined by the frozen 
reservoir. Analytical evaluation of the trace within Eq. <0 
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leads to the following expression for the nonequilibrium pop¬ 
ulation dynamics, 

/ g~ + A 2 cost t't ) 

dVd(± p bath (P, Q, 0)- g2 + A , , (15) 

where f = 2 Vs 2 + A 2 is the Rabi frequency of the modi¬ 
fied system Hamiltonian. The integration over the ensemble 
of equilibrium configurations of the bath reduces to averag¬ 
ing over different values of T c, (0) that are consistent with the 
bath distribution function. For some realizations of the bath, 
Eq. 0 recovers the Rabi oscillations characteristic of the 
isolated system if A cl ( 0) = 0. Conversely, when T c, (0) + 0, the 
population starts from 1 at t = 0 and will oscillate between 
(e 2 + A 2 )/(§ 2 + A 2 ) = 1 and (e 2 - A 2 )/(e 2 + A 2 ). Taking for 
simplicity s — 0, one notes that the lower bound of the popula¬ 
tion oscillations increases with increasing A cl ( 0), approaching 
1 as A c, (0) —> oo. This limit corresponds to an infinitely rigid 
bath that completely localizes the excitation on its initial site. 

Averaging over different realizations of the slow modes de¬ 
creases the amplitude of oscillations in the population dynam¬ 
ics due to the decoherence between the functions with distinct 
oscillation frequencies. In general, Redfield theory has dif¬ 
ficulty describing non-Markovian, multi-step relaxation dy¬ 
namics. However, when T s i ow = iC x f Q doj J s i ovl (u))/u> is suf¬ 
ficiently large, the full dynamics produced by Redfield theory 
with frozen modes at finite u>* will include both a slow, per¬ 
haps oscillatory component as well as a more rapid decay in¬ 
duced by the high frequency modes in Jf ast (to). These qualita¬ 
tive considerations suggest this approach may correct certain 
deficiencies of conventional Redfield-like approaches. In the 
next section, we test the approach quantitatively. 


III. RESULTS 

In the following, we compare the numerically exact popu¬ 
lation dynamics reported by Thoss et alW^ for the spin-boson 
model with a Debye spectral density and the initial condi¬ 
tion T(0) = |1) (1| exp(-y8// batb )/Z batb with the results obtained 
from the Redfield-FM method. Subsequently, we examine 
the effect of relaxing the mode-freezing approximation by 
treating the low-frequency modes dynamically via the Ehren- 
fest method. We call this latter approach the hybrid-Redfield 
method, in analogy with the previously developed hybrid- 
NIBA method.® 


A. Computational Details 


To treat the frozen portion of the spectral density, J s \ ow (tS), 
we have discretized the bath in to f = 300 modes with fre¬ 
quencies and couplings given bySLUJ 


and 


u>k = u) c tan 


2/ ( *- 1/2) 


2A 


f 


(16) 


(17) 


(a) aj c = 0.25, A = 0.25 ,/3 = 5.0 (b) oj c = 0.25, A = 0.25, {3 = 0.5 




Time t [A 1 ] Time t [A *] 


FIG. 2. Results from Redfield-FM approach compared with standard time- 
local Redfield theory and exact numerics, of = max[o) c , ^ | and cjr = 
2 VA- + e 2 . All units are scaled by the electronic coupling A. All panels 
correspond to unbiased cases (s = 0), except for panel (b), where e = 1.0. 
Other parameters are stated in the panels. 


Initial conditions for the reservoir of frozen modes were 
sampled from a Wigner distribution. Sampling from this dis¬ 
tribution becomes particularly important at very low tempera¬ 
tures, where quantum effects become significant. However, 
for most cases, sampling from a Boltzmann distribution is 
sufficient since the modes being samples are always low- 
frequency. For convergence, up to 10 4 trajectories have been 
run for the results presented. 


B. Redfield-FM Method 


As mentioned in Sec. IIB the validity of Redfield theory 
is limited to the small r\ regime. Fig. 2(a) shows the results 
for a slow bath ( a> c = 0.25) with small reorganization en¬ 
ergy (/l = 0.25) at low temperature (J5 — 5.0); here and in 
the following, all energies are in units of A. In spite of the 
slow bath, the dimensionless applicability parameter is only 
slightly larger than unity (tj = 1.6) suggesting that Redfield 
theory should be reasonably accurate, in agreement with the 
numerical results. The Redfield-FM method provides an even 
better estimate of the dynamics, almost quantitatively correct¬ 
ing the already accurate Redfield dynamics. 

Fig. 2(b) considers a biased system (s = 1.0) with the same 
parameters, except at much higher temperature {J3 = 0.5), 
yielding an applicability parameter which is now significantly 
larger than unity (rj = 16). Here, it is evident that the Redfield 
dynamics relax far too quickly, suppressing the coherence and 
missing the slower relaxation process revealed by the exact 
dynamics. The improvement afforded by the Redfield-FM 






















5 


method compared to standard Redfield theory is clear. In par¬ 
ticular, the Redfield-FM approach accurately reproduces the 
short- to intermediate-time dynamics, the frequency of the 
oscillations, and the initial rate of decoherence. The terminal 
decay rate is slightly underestimated due to the mode-freezing 
approximation, as discussed in Section II.B. Yet despite these 
shortcomings, the improvement derived from a simple scheme 
like the Redfield-FM approach with the numerical complexity 
of the original Redfield theory is noteworthy. 

For cases exemplified by Fig. 2(c), serious problems such 
as the violation of the positivity of the RDM dynamics can oc¬ 
cur within standard (non-secular) Redfield theory. Figure 2(c) 
corresponds to a slow bath (oj c = 0.25) and a large reorgani¬ 
zation energy (A = 5.0) again at high temperature (J3 = 0.5), 
for which the applicability parameter is very large (77 = 320). 
Despite the evident failure of the Redfield equations to even 
maintain positivity, the Redfield-FM method is able to cor¬ 
rect the positivity issue and almost quantitatively reproduce 
the two-step relaxation process in the exact dynamics up to 
intermediate times. 

It is possible to understand the surprising success presented 
in Fig. 2(c) in the context of the analysis of Section II.B. Using 
the definitions given in that section, the effective parameters 
for the Redfield equation are Tf ast = 2.5, a> c = 0.5, yielding 
77 = 40. Although 77 » 1, the reduction by an order of mag¬ 
nitude from the initial value, 77 = 320, is sizeable, and likely 
responsible for solving the positivity problem evident in the 
bare Redfield dynamics. The reproduction of the two-step re¬ 
laxation process is a direct result of the trapping effect that 
arises from freezing a large portion of the low-frequency bath 
in the presence of large coupling. This example indicates that 
the trapping effect can partially reproduce slow relaxation dy¬ 
namics associated with strong system-bath interactions. 

Figure 2(d) shows the regime of intermediate bath speeds 
(cj c = 1), large reorganization energy (A = 2.5), and inter¬ 
mediate temperature (J3 — 1). In contrast to Fig. 1(c), the 
Redfield-FM method is not capable of significantly improv¬ 
ing the Redfield dynamics in this regime, missing the two- 
step relaxation process visible in the exact dynamics. In light 
of the previous case, it is evident that the slowing down of the 
RDM dynamics can be caused by freezing a large portion of 
the strongly coupled modes, an effect which is absent in this 
case. In this case, /lf ast = 2.1 and T s i ow = 0.4, which indicates 
that most of the reorganization energy is included already in 
the high-frequency portion of the bath. In such cases, the 
Redfield-FM method will yield results that are similar to bare 
Redfield theory. 

We now address the dependence of the dynamics on the 
choice of of . Eschewing the simple criteria for choosing of 
presented above, one may ask how closely the Redfield-FM 
dynamics can be made to agree with exact dynamics when of 
is allowed to vary. To address this question, we include two 
extreme cases in Fig. 3. First, Fig. 3(a), which corresponds 
to the same parameters as those of Fig. 2(c), shows that opti¬ 
mization of of can result in quantitative agreement between 
the Redfield-FM result and the exact dynamics. Such agree¬ 
ment may be understood as the result of fortuitous cooperation 
between strongly dissipative Redfield dynamics that damp the 


(a) w c = 0.25, ,1 = 5.0, y3 = 0.5 (b) ai c = 1.0, A = 2.5,/3= 1.0 



FIG. 3. Redfield-FM results for tu* = max[tu c , ^ ] and an “optimized” value 
for tif. Both cases considered here correspond to e = 0 and all units are scaled 
by the electronic coupling, A. Note that the set of parameters for panels (a) 
and (b) in this figure correspond to the set of parameters in Fig. 2, panels (c) 
and (d), respectively. 


frozen mode-generated oscillations and the trapping effect 
from the mode-freezing approximation that prevents imme¬ 
diate relaxation to the equilibrium population. Conversely, 
Fig. 3(b), which corresponds to the parameters in Fig. 2(d), is 
an example of when perfect agreement is impossible. Clearly, 
attempts at optimizing of result in better agreement of the 
two-step relaxation process at the cost of long-lived oscilla¬ 
tions, a direct result of including a large fraction of modes into 
the slow part of the bath. In freezing a sufficiently large por¬ 
tion of the reservoir to reproduce the trapping effect, /lf ast is re¬ 
duced to the point where the Redfield dynamics are no longer 
sufficiently dissipative to damp the frozen mode-generated 
oscillations. In addition, /l s i ow = n J () do>J s \ ow (oj)/oj is also 
not large enough to ensure that the oscillations dephase suf¬ 
ficiently rapidly. Overall, it is clear that although it may be 
possible to optimize the results, the simple initial criteria pre¬ 
sented represent a robust approach to frozen mode dynamics 
that essentially always yields results that are as good or better 
than bare Redfield dynamics without a significantly increased 
computational cost. 


C. Relaxing the Mode-Freezing Approximation: A 
Dynamical Hybrid Redfield Method 

On first inspection, the mode-freezing approximation ap¬ 
pears extreme. To thoroughly assess its effect, we develop a 
dynamical hybrid method in which we evolve the previously 
frozen low-frequency modes in J s i ov/ (of) via classical Ehren- 
fest dynamics. The derivation and the implementation details 
of this approach may be found in Appendix |C| 

This hybrid-Redfield method is similar in spirit to the suc¬ 
cessful hybrid-NIBA developed and implemented in Ref. |9] 
Evolution of the low-frequency modes using Ehrenfest dy¬ 
namics in such hybrid approaches only requires two assump¬ 
tions: (i) that r(0 » Psiow(0Psys+fast( ? )> and (ii) that the motion 
of the low-frequency modes is well-captured by classical me¬ 
chanics. For such a factorization to be valid, the reorganiza¬ 
tion energy due to the low-frequency bath needs to be small, 
i.e., /1 S ] 0W = 7 r~ 1 f () doj J s \ ow (ui)/oj < A. The applicability 
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of classical dynamics relies on the low energies of the reser¬ 
voir modes and sufficiently high temperatures that help sup¬ 
press quantum effects.® 2 However, even when the Ehrenfest 
approximation is valid, problems may arise. Most prominent 
among these is that the final populations approach those of the 
infinite-temperature limit. 12 

In contrast to the Hamiltonian derived under the mode¬ 


freezing approximation in Eq. (13i, the modified Hamilto¬ 
nian that needs to be treated via the Redfield equation in the 
hybrid-Redfield method is time-dependent, 


■ Act,. 


H"(t) = [e + A d (t)}cr z 

+ O'- 2 CkQk + \ Z + W *G*]> 


(18) 


ke fast 


where the disorder due to the low frequency bath is no longer 
static as it is in the Redfield-FM method, but rather dynamic, 
namely A c, (t) = £ A . c k Q k (t). 

Since the system part of this Hamiltonian is nondiagonal 
and time-dependent, evolution with respect to the system 
Hamiltonian requires diagonalization at every time-step, sig¬ 
nificantly increasing the computational cost associated with 
the method proposed here. The need to evolve the low- 
frequency bath also adds to the computational cost of the ap¬ 
proach. Importantly, under the mode-freezing approximation, 
we circumvent these costly requirements. This means that, 
aside from the trivial cost of parallelization for the ensemble 
averaging over the slow bath, the Redfield-FM method scales 
as gracefully with system size as the original Redfield equa¬ 
tion. 

For completeness, we remark that the nonequilibrium pop¬ 
ulation dynamics under the hybrid-Redfield approximation 
now take the form 

(0"-(O) ~ I dPdQ p s iow(P, Q, 0 Tr sys+fast [p sys+fast (f)cr,]. 

(19) 

Fig. 4 shows two sets of parameters for which the hybrid- 
Redfield scheme yields results that illustrate the issues at play 
in comparing the hybrid-Redfield approach to the Redfield- 
FM method. Extensive testing of the hybrid method suggests 
that an approximately optimal form for the splitting frequency 
can be taken as 


- wrT/w,.. (20) 

Physically, this form encodes the interplay between the Red¬ 
field and Ehrenfest methods, favoring a larger portion of the 
modes to be treated classically with increasing Rabi frequency 
u>r, which is a measure of how rapidly the electronic system 
evolves. Furthermore, this form of of ensures that in the limit 
of small A and large u> c , the hybrid method correctly repro¬ 
duces the more appropriate Redfield dynamics, whereas in the 
limit of large A and small a> c , it reproduces the Ehrenfest re¬ 
sults. It is expected that generally nontrivial results may be 
obtained from this method for cases where of ~ u> c , as is the 
case for the choice of of used in the Redfield-FM approach. 

Fig. 4(a) corresponds to the parameters in Fig. 2(c) and il¬ 
lustrates that, by means of the suggested form for of, hybrid- 


(a) (D c = 0.25, ,1 = 5.0, y3 = 0.5 (b) ai c = 1.0, 4 = 2.5 ,j8= 1.0 



Time t [A 1 ] Time t [A 1 ] 


FIG. 4. Hybrid-Redfield results for co*^ = cjr . Both cases considered here 
correspond to e = 0 and all units are scaled by the electronic coupling, A. 
Similar to Fig. 3, the set of parameters for panels (a) and (b) in this figure 
correspond to the set of parameters in Fig. 2, panels (c) and (d), respectively. 


Redfield automatically tunes itself to yield nearly optimal re¬ 
sults achievable from the two parent methods. This exam¬ 
ple, for which of h » o> c , illustrates that the hybrid-Redfield 
method trivially reproduces the Ehrenfest result when it is ap¬ 
propriate. It is noteworthy that the Redfield-FM method ob¬ 
tains similar agreement at a much lower computational cost 
without evolving the reservoir modes, indicating that dynamic 
treatment of these modes may not be generally necessary. In¬ 
deed, it is rather remarkable that the Redfield-FM approach 
basically recapitulates the Ehrenfest results even though no 
Ehrenfest dynamics are used. 

Fig. 4(b) shows the analogue of Fig. 2(d), where the 
Redfield-FM method fails to correct the Redfield dynamics. 
In contrast, the hybrid-Redfield results are in very good agree¬ 
ment with the exact dynamics. Indeed, the hybrid method is 
able to qualitatively and almost quantitatively reproduce the 
shape of the two-step relaxation process evident in the exact 
dynamics, an effect that both Ehrenfest and Redfield dynamics 
independently miss. 

As the above considerations indicate, there are cases where 
the dynamical hybrid-Redifeld method can provide a substan¬ 
tial improvement over the Redfield-FM method, albeit at a 
much higher computational cost. In most regions of parameter 
space we have studied, however, we find that hybrid-Redfield 
theory offers little accuracy gain over the Redfield-FM ap¬ 
proach. Thus, the benefits of the hybrid-Redfield approach 
compared to the Redfield-FM method do not justify its use 
when accuracy and cost are factored together. 


IV. CONCLUSIONS 

In this work, we have presented a new scheme for simulat¬ 
ing dynamics in quantum dissipative systems. Our approach, 
which we call the Redfield-FM method, recognizes that stan¬ 
dard Redfield theory becomes inaccurate for slow bath de¬ 
grees of freedom. By partitioning the bath into high- and 
low-frequency components, we propose solving the Redfield 
equations for the high-frequency partition in the statically 
disordered field of the low-frequency components. Such an 
approach may greatly increase the accuracy of Redfield the- 
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ory in highly non-Markovian regimes at essentially the same 
computational cost. In addition, we find that this simple ap¬ 
proach can fundamentally cure positivity problems associated 
with standard non-secular Redfield theory. We have further 
discussed a scheme (the hybrid-Redfield approach) whereby 
the previously frozen degrees of freedom are instead evolved 
with classical Ehrenfest dynamics. While this method can im¬ 
prove upon the dynamics as described by the Redfield-FM 
approach, the increase in accuracy is incremental and comes 
at a significantly larger computational cost. Overall, while the 
Redfield-FM method does not cure all of the ills of Redfield 
theory, it does provide a simple and efficient framework for 
improving its accuracy and range of validity, especially for 
sluggish bath degrees of freedom such as those implicated in 
biological energy transfer. 
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Appendix A: Derivation of Redfield equations 


following exact equations of motion, 

ymo = -iPZiitW + Q)Y(t) (A2) 

at 

JrQYii) = -iQ£/(t)(P + Q)Y(t). (A3) 

at 

Formal integration of Eq. ( |A3[ ) yields 

Qrj(t) = -i f dr g(t,T)QJL I (typr I (r), (A4) 
Jo 

X f 

ds <3i)/(s)] and the time order¬ 
ing (+) implies that time arguments increase from right to 
left. Substitution of this expression in Eq. ( |A1[ ) results in the 
Nakajima-Zwanzig equation,EMI which is expressed in terms 
of the time convolution of a memory term with pi(t) at earlier 
times as 


j-P/M = - f dT Kit - T)p/(r), (A5) 

dt Jo 


where K(t - r) = Tr balh [£i(t)g(t, T)<3£,(T)p ba th] is the (time- 
nonlocal) memory function. 

If, instead, we use the formal solution of Eq. (A1 1 to evolve 
E Jt) backwards in time to an earlier time r, we obtain 


Here, for completeness, we review the derivation of the 
Redfield equations. For a more detailed discussion of Red¬ 
field theory, we refer the reader to Refs. fl3]and fl4[ 

In the following development we utilize a projection op¬ 
erator technique to derive an equation of motion for the re¬ 
duced density matrix (RDM) of the system, defined as pit ) = 
Tr bat h[T(f)], where T(f) = e^' H, YiO)e' Hl and HO) is the ini¬ 
tial density matrix of the full system and bath. Moreover, we 
assume that the initial condition for the (total) density ma¬ 
trix contains no system-bath correlation, such that T(0) = 
p(0)pbath- where p(0) is an arbitrary Hermitian system oper¬ 
ator, p bath = e^ Hbath /Z, Z = Tr b ath[e“ /3Wba,h ] and p = 1 /k B T is 
the thermal energy. Treatment of general initial conditions is 
also possible via the projection operator technique at the ex¬ 
pense of the introduction of additional inhomogeneous terms 
in Eq. (j^EE [ n the following, we ignore initial correla¬ 
tions, but note that their inclusion in the present framework 
is straightforward. 

We start from the Liouville equation for the full density 
matrix in the interaction picture where the total Hamiltonian 
is divided into a zeroth order part and an interaction part, 
H — H 0 + Hi, such that 

A/(0 = -iZMYJt), (Al) 

dt 

r 7 (0 = e iH °T(t)e~ iHot and £ 7 (f) = [e~ iHQt H x e iH ^\ ...]. To ob- 
tain the dynamics of the RDM, we define a projection oper¬ 
ator of form P = Pbath Trbath [ • ■ • ] with Q = 1 - P. We note 
that action of P on the full density matrix followed by trace 
over the bath results in the RDM in the interaction picture, 
pi(t) = TrbathlTT/O)]. Using these definitions, we obtain the 


T/(r) = Git^YJt), 


(A 6 ) 


X f 

ds£i{s)\ and the time ordering (-) 
requires that time arguments increase from left to right. Re¬ 
placing this expression in Eq. (A4i and solving for <2E/(f) 
yields 


QYjt) = [i - mr'mPYjt), (A7) 


where 


-■/' 


= -i dr git, T)Q£(T)PG(t, t). 


(A 8 ) 


We note that a crucial requirement for the validity of this 
derivation is the existence of [1 - 2 (/)] -1 . 

Substitution of Eq. (A7 1 into Eq. ( |A1| > and subsequent trace 
over the bath degrees of freedom results in the following time- 
local equation of motion for the RDM,El 


-rPiit) = Kit)pit), 
dt 


(A9) 


where R(t) = —f Tr ba th[X/(T)[ 1 - 1(f)] 'Pbath] is the (time- 
local) rate function. 

The expression for the dynamical evolution in either the 
time-nonlocal (Eq. ( |A5| >) or time-local (Eq. ( |A9| >) form is ex¬ 
act but prohibitively difficult to evaluate without resorting to 
approximation schemes, such as truncated generalized cumu- 
lant expansions. Perturbative expansion to second order in the 
system-bath coupling (where H\ -V from Eq. ([3])) results in 
a non-Markovian generalization of the Redfield theory. Alter¬ 
natively, one may derive both forms of the Redfield equations 
via resummations of differently time-ordered cumulants.E 19 ! 
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These derivations explicitly show that both forms of gener¬ 
alized Redfield theory account for non-Markovian behavior 
and have similar applicability require mentsPS 21 Specifically, 
since Redfield theory is tantamount to second-order pertur¬ 
bation theory in the system-bath coupling, truncation at low 
order is only accurate for rj < 1 , where // = max[Jj|, ■££-] 
is the validity parameter introduced in Sec. IIB Despite this 
restriction, the Redfield equations have been shown to per¬ 
form surprisingly well, often beyond the small-T and large- 
oj c regimesP 223 Nevertheless, for inappropriate regions of pa¬ 
rameter space, severe problems can arise, such as violation of 
positivity in the reduced density matrix . 3 


Appendix B: Markovian Redfield theory with Frozen Modes 

We wish to consider the performance of the frozen modes 
method for a strictly Markovian version of Redfield theory, i.e. 
with a rate tensor R = R(t —» oo). In this limit, the time inte¬ 
grals become Fourier-Laplace transforms, such that the Red- 
field tensor elements can be expressed algebraically in terms 
of the spectral density J(a>) evaluated at energy differences 
hoJij = (Ej - fjjP 4 23 More specifically, we are interested in 
the dephasing terms of the Redfield tensor, which in general 
contain an elastic contribution 

Rijij ~ 8 2 ,j = 0 +) hbe(w = 0 +). (Bl) 

At low frequencies, the Bose-Einstein distribution, nm-ico) ~ 
kT/u>, such that for a spectral density of the form J(u>) ~ of, 
we find 

R liu ~ kTto s ~ l \ . (B2) 

11 L=o 

For all ‘super-Ohmic’ spectral densities with s > 1, this elas¬ 
tic contribution to the dephasing rate vanishes. However, for 
an Ohmic spectral density with s = 1 there is a pure dephasing 
rate which vanishes only at T = 0. This contribution to the de¬ 
phasing rate in the system’s eigenbasis can significantly affect 
both the population and coherence dynamics in the original 
basis of the problem. 

We now return to the idea of a frozen modes variant of 
Markovian Redfield theory. Consider specifically an Ohmic 
spectral density with any non-zero splitting frequency of. 
After partitioning, the fast spectral density has the low- 
frequency behavior 7f as t(w) ~ of with s > 1 , which yields 
no elastic contribution to the dephasing rate. For this reason, 
a frozen modes version of Markovian Redfield theory does not 
reduce to the Redfield limit until the singular point of = 0. In¬ 
stead, as of —> 0, the result approach a Redfield result which 
neglects the Ohmic pure dephasing rate. We emphasize that 
the time-dependent variants of Redfield theory are not sign- 
ficantly affected by this problem until very long times, and 
that all methods are only affected for strictly Ohmic spectral 
densities. 

We propose a very simple solution to this pathological be¬ 
havior in Markovian Redfield theory, by modifying the fast 
spectral density via 

■/ fast (w) = [1 - S(a>, u>*)]J(co) + W(oj, e)J(a>), (B3) 



(c) A = J,r?= 15 (d) A = 5J, t] = 75 



Time t [fs] Time t [fs] 


FIG. 5. Comparison of numerically exact (HEOM) population dynamics to 
the results of standard Markovian Redfield theory, a straightforward variant of 
Markovian Redfield theory with frozen modes (Red-FM), and a dephasing- 
corrected variant as discussed in the text (Red-FM-D). The system-bath 
Hamiltonian is that of Ref.j4]with e = 100 cnT 1 , J = 100 cnT 1 , = 100 
fs, and T = 300 K. 


where W(o>, e) is a rectangular window function centered at 
the origin with width e, and e should be chosen very small. 
In this way, the ‘fast’ part of the bath will always produce a 
pure dephasing rate for arbitrary splitting frequency of. Thus, 
the Markovian Redfield-FM dynamics will smoothly inter¬ 
polate towards the standard Markovian Redfield result as cf 
approaches zero. In Fig. 5, we compare the results of stan¬ 
dard Markovian Redfield, Markovian Redfield-FM without 
this dephasing correction, and Markovian Redfield-FM with 
the correction. Results are presented for the model excitonic 
dimer discussed by Ishizaki and Fleming.® Importantly, we 
find that this correction typically improves the results of the 
Markovian Redfield-FM variant, quite significantly in cases 
of strong system-bath coupling. 


Appendix C: RDM Hybrid Method 

Here we relax the mode-freezing approximation by deriv¬ 
ing a fully hybrid method that separates the complete system 
into a slow part consisting of the low-frequency component 
of the bath, and a rapidly-evolving part that includes both 
the electronic system and the high-frequency portion of the 
phonon bath. In this hybrid scheme, the slow part is treated 
quasi-classically, while the fast part is treated at the level of 
Redfield theory. Overall, the fast (slow) component of the 
system evolves in the mean field of the slow (fast) one. 

Other hybrid approaches that combine classical and quan¬ 
tum dynamics include the self-consistent hybrid method of 
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Wang and coworkers®2D which yields numerically exact 
dynamics, and the approximate hybrid-NIBA approach of 
Refs. [8jand[9]In the former method, at*, which is the energy 
scale that determines the splitting of the bath into slow and 
fast parts, is strictly a convergence parameter. In the latter, of 
is an empirically determined adjustable parameter. As an ap¬ 
proximate method, the hybrid-Redfield scheme derived here 
is akin to the hybrid-NIBA method. For a more detailed dis¬ 
cussion of the hybrid RDM method, we refer the reader to 
Ref. E 

As in the Redfield-FM approach, we make the approxima¬ 
tion that r(0 as PsiowPsys+fastM. where p sys+fast (f) is the den¬ 
sity matrix for the system and fast bath degrees of freedom 
and PsiowW is the density matrix for the slow bath degrees of 
freedom. The system and fast bath modes obey the following 
effective Liouville equation, 


d(psys+fast(0 

dt 


i[H (0>Psys+fast(0L 


(Cl) 


where 


bath correlation, as in the case of the Redfield-FM method, 
takes the following form, 

1 f°° 

C(t) = - du> 7fast(w) 

n Jo (C6) 

x [coth(j6w/2) cos(cuf) - /sin(wf)]- 


To obtain the results shown in Fig. 3, trajectories cor¬ 
responding to a set of initial conditions sampled from the 
Wigner distribution 27 were calculated via a second-order 
Runge-Kutta scheme, using a step size of dt = O.OIA -1 . As 
required by the Runge-Kutta procedure, if ft) was kept con¬ 
stant during the evolution of the bath while A cl (t) was kept 
constant during the evolution of the system. Explicitly, over a 
half-time step, the equations for the bath become 



(C7) 


H"(t) = [s + A c (t)]cr z + Acr x 


+ CkQk + \ Yj \ p2 k + ^2*]’ 


(C2) 


ke fast 


and 


Pk[t+ y) = Pk(t) cos+ u k a k (i) s in, ( C8 ) 


and A cl is a dynamically fluctuating bias, A cl (t) = 

Xfcslow CkQk(t')- 

A classical treatment of the reservoir leads to the following 
equations of motion. 


and 


dQk 

dt 


= Pk, 


(C3) 


^ = -<J k Qk ~ c k d- z (t). (C4) 


Employing the Ehrenfest approach demands that each part of 
the system evolves in the mean field of the other. For the 
quantum portion, the classical mean field consists of the time- 
dependent contribution to the bias, A cl (t). For the classical 
portion, the force term -c k & z {t) = -c k ^^+^[^^,.+(^(0] 
in the equations of motion embodies the mean-field ‘back- 
reaction.’ This force term moves the classical oscillators from 
the ground state minima to the displaced minima associated 
with the excited state. 

Using Eqs. (|5]» and ( |C2[ >, the time-local Redfield equation 
takes the form 

j t p(t) = -i[Hj y ,(t),p(t)] 

ds {C(s)[cr z ,cr-(-?)p(f)] ( C5 ) 

- C*(s)[cr z ,p(t)cr : (—i)]}. 



where cr z (t) = U ] Q (t)o- z U 0 (t ), U Q (t) = exp[-i ^dr W" s (t)], 
and = [s+A c \t)]cr z + Acr x + Zkesio»[P 2 k+u 2 k Ql]/2. The 


where 

a k (t) - Q k (t) + ^<r z (t). (C9) 

In the hybrid-NIBA method of Ref. [9] the zeroth-order 
propagator necessary to evolve the perturbation in the interac- 
tion picture, Uif t) = exp[-/ J c/rf/" s (T)], was simple to cal¬ 
culate since //" s was diagonal. In contrast, /-/'' s (f) for hybrid- 
Redfield contains off-diagonal elements. Within the Runge- 
Kutta scheme, this obstacle is easy to overcome, though it re¬ 
quires diagonalization of the time dependent H' S y S (t) at every 
time step. Because numerical diagonalization at every time- 
step is necessary for systems with more than two degrees of 
freedom, this can become computationally expensive for suf¬ 
ficiently large systems. 
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